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Quantum vortex dynamics in Bose-Einstein condensates or superfluid helium can be informa- 
tively described by the Gross-Pitaevskii (GP) equation. Various approximate analytical formulae 
for a single stationary vortex are recalled and their shortcomings demonstrated. Significantly more 
accurate two-point [2/2] and [3/3] Pade approximants for stationary vortex profiles are presented. 
Two straight, singly quantized, antiparallel vortices, located at a distance do apart, form a vortex 
dipole, which, in the GP model, can either annihilate or propagate indefinitely as a 'solitary wave'. 
We show, through calculations performed in a periodic domain, that the details and types of be- 
havior displayed by vortex dipoles depend strongly on the initial conditions rather than only on 
the separation distance do (as has been previously claimed). It is found, indeed, that the choice of 
the initial two- vortex profile (i.e., the modulus of the 'effective wave function'), strongly affects the 
vortex trajectories and the time scale of the process: annihilation proceeds more rapidly when low- 
energy (or 'relaxed') initial profiles are imposed. The initial 'circular' phase distribution contours, 
customarily obtained by multiplying an effective wave function for each individual vortex, can be 
generalized to explicit elliptical forms specified by two parameters; then by 'tuning' the elliptical 
shape at fixed do, a sharp transition between solitary- wave propagation and annihilation is captured. 
Thereby, a 'phase diagram' for this 'AnSol' transition is constructed in the space of ellipticity and 
separation and various limiting forms of the boundary are discussed. 



I. INTRODUCTION 

The Gross-Pitaevskii (GP) equation^"— is widely ac- 
cepted basic model to study vortex dynamics in 
superfluid helium or Bose Einstein condensates. It has 
the merit of providing an effective quantum mechanical 
description of the vortex core dynamics, while also al- 
lowing detailed predictions of vortex reconnection^"— a 
phenomenon which is thought to centrally affect the large 
scale behavior of quantum turbulencci ^^i^^ 

The general focus of this present work is the specifica- 
tion of appropriate initial conditions to perform system- 
atic, reproducible, vortex dynamics calculations through 
the GP equation. This is a matter of some delicacy, be- 
cause given the nonlinearities involved one must rely al- 
most entirely on numerical solutions, which, as we will 
demonstrate, depend strongly on the details of the ini- 
tially imposed complex, efFcctive wavcfunction ^'(x, t — 
0). Of course, the boundary conditions also play a role, 
although for the present purposes, a less crucial one. In 
this study we will always assume periodic boundary con- 
ditions. 

We address the general objective here by study- 
ing single vortex profiles and vortex dipole dynamics, 
which may be usefully relevant to physical situations 
that can be regarded as two-dimensional, such as he- 
lium filmSfi^ or experiments performed in trapped Bose- 
Einstein condensates!^"— 

In order to select initial conditions that lead to simu- 
lated vortex dynamics, two conceptual steps are needed: 
first, one needs to describe mathematically a quantum 
vortex appropriately for the GP equation; second, one 



needs to set up one or more vortices in the computational 
domain. 

This second step has usually been achieved merely by 
multiplying together wave functions describing each in- 
dividual vortex4i^iiii^ For exploratory purposes this de- 
vice is simple and may well be adequate. However, if one 
wishes to address more quantitative and subtle issues, 
such as arise naturally in contemplating the wealth of 
experimental dat a^^i^^'^° and what may be found in the 
future, a more systematic approach is called for. Here 
we describe a different attack which exploits the prop- 
erties of the diffusion equation associated with the GP 
equation and thus allows one to impose low-energy ini- 
tial conditions. The use of the diffusive GP equation is 
sometimes referred to as the imaginary time propagation 
method^^^^ 

The mathematical description even of a single vortex 
entails some complications. A single straight vortex in 
an infinite domain is, in fact, optimally represented as a 
solution of the stationary GP equation. However, there is 
no convenient analytical form to express it precisely; as a 
result, a numerical study is required. However, this may 
not be convenient or practical when wave functions for 
multiple vortices are required. Thus it has been custom- 
ary to employ some analytically convenient but intrinsi- 
cally approximate profiles. Various approximate profiles 
for a single straight vortex have been proposed. We recall 
three of these, namely, in order of increasing complexity, 
the Fetter;^ Kerr— and Berlof^ approximants; how- 
ever, as distinct improvements we propose new [2/2] and 
[3/3] Pade approximants. It is worth stressing, in this 
regard, that the GP equation conserves energy. Con- 
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scqucntly any energy imposed by means of vortex ini- 
tial conditions that exceeds the minimal energy entailed 
in the exact vortex configuration is not dissipated but, 
rather, dispersed by mechanisms which are not easy to 
predict a priori, and often confusing to interpret. 

In the GP equation unit vortex dipoles, characterized 
by a separation distance do, can either annihilate and 
emit their energy via outgoing waves, or, after a tran- 
sient, propagate at constant velocity in the form of soli- 
tary waves. Furthermore, Ivonin^^ and Ogawa et al^ 
claimed there was a nonzero critical value of do, be- 
low which no solutions describing dipoles moving uni- 
formly were possible. However, Jones and Roberts^^ and 
BerlofS^i^ found accurate numerical solutions for con- 
stant velocity solitary waves at separations down to the 
limit^ do = 0, which corresponded to a propagation ve- 
locity U — 0.45, in dimensionless units where the speed 
of sound iai^ c = l/-\/2- As we show, this discrepancy 
highlights the importance of the initial vortex configura- 
tions in studying solutions of the GP equation. In partic- 
ular, the phase pattern of the initial order parameter is 
significant. Indeed, by varying it, we are able to generate 
solutions in which the vortex dipoles either annihilate or 
propagate. Our GP calculations also serve to check the 
dependence of the time-to-annihilation and of the dipole 
propagation velocity, U, on the initial separation as do 
^ 0. 

In Sec. II we recall the Gross-Pitaevskii and we de- 
scribe the Diffusive Gross-Pitaevskii (DGP) equation. 
Section III is devoted to the study of single vortex pro- 
files, while in Sec. IV we present some illustrative vor- 
tex dipole calculations that demonstrate how different 
choices of even very similar initial conditions significantly 
affect aspects of the observed phenomenon. Finally our 
results are summarized in Sec. V. 



II. THE GROSS-PITAEVSKII EQUATION 

The Gross-Pitaevskii equation,—"— expressed in terms 
of an effective complex bosonic or 'condensate' wave func- 
tion ^'(x, t), where x is a d{=2 or 3)-dimensional spatial 
coordinate with t denoting the time, is 

^h^^-—W^^ + Vo^\^\'-^^^. (1) 

Here 27Th is Planck's constant, m is the mass of the Bose 
particles being simulated, while /i is the chemical po- 
tential and Vo represents the strength of the short-range 
boson-boson repulsive potential. Evidently, the nonlin- 
ear term accounts for the interatomic interactions. 

This equation can be made dimensionless by rescal- 
ing the coordinates x, t, and ^'(x,i), respectively, by a 
characteristic length, the healing length^ = h/y/2m^ 
(of order 0.5 A for helium-4), by a characteristic time, 
to ~ m^'^/h (of order 10~^ ps for helium-4), and by a 
characteristic modulus, chosen as the value of the wave 
function in an infinite domain prior to the onset of a 



disturbance^^ l^ool = \/ m/M]- The dimensionless equa- 
tion then becomes 

-2i^ = V2* + (l-|*|2)*. (2) 

The total energy, E, associated with this form, measured 
inside a domain VL with respect to the uniform state 
^'oo = 1 is conserved, and is given by the sum of the 
kinetic energy Ek, and the potential energy Ei with^^ 

EK = lJjy^^d^ (3) 

and 

Ej^^J^il~\^\Yd^. (4) 

The kinetic energy can be further split into a classical 
kinetic energy, Eck, and a quantum energy, Eq, which 
are given by, respectively^^ 

EcK^\j f\V^?d^ (5) 

where is expressed as ^ = /(x, <)e"^*^'''*\ and 

EQ^\j^{V.ffd^. (6) 

It is useful, for thinking in physical terms, to define the 
density of the condensate p as p = f^, and the velocity 
field V as V = V(/). These equivalences allow to write the 
GP equation in a hydrodynamical form.— 

The diffusion equation associated with the GP equa- 
tion - the DGP equation - is 

2^ = V2* + (I- (7) 
at 

This does not conserve energy; rather, when t ^ oo, one 
has 'J — > ^, where ^ minimizes the total energy and 
coincides with a stationary solution of the GP equation. 
Therefore, the DGP equation can be used to find fixed 
points for the GP equatioui^ (As noted above, this imag- 
inary time propagation method, has been often used by 
the BEG community. 

However, even if the focus is rather on the dynamical 
behavior of unstable structures, such as vortex dipoles or 
antiparallel vortices with three-dimensional geometries, 
the DGP equation can be conveniently exploited to find 
low-energy initial conditions. Here this is demonstrated 
for vortex dipoles. 

In our calculations)^ the GP equation has been 
integrated numerically through a split-step spectral 
method^ in a periodic domain. The use of the Dis- 
crete Gosine Transform (DGT), rather than the Dis- 
crete Fourier Transform (DFT), allows one to account 
for periodic images by calculations in a domain smaller 



3 



by a factor of four (in two dimensions). The DCT, in- 
deed, implies an even extension of the original function 
which saves us from having to add three reflections of 
each single vortex in the computational box, a condi- 
tion otherwise essential to meet the appropriate periodic- 
ity conditions.— The initial conditions (before any relax- 
ation with the DGP equation) are obtained by multiply- 
ing the wave-function that represents the vortices inside 
the domain by the wave-functions that correspond to the 
first eight periodic images outside the domain to ensure 
that ^ at the boundaries is sufficiently smoothi^ii^i^ 



III. VORTEX PROFILES 

In the language of the GP equation, a vortex line is a 
long, thin, locally cylindrical structure, straight, bent or 
wavy, defined by a linear locus of zeros of the effective 
wavefunction. In particular, a quiescent straight vortex 
is represented by a stationary axisymmctric solution of 
(PI which can be written in cylindrical coordinates, (r. 
(j), z), as 



^{r,(f>,z) = f{r) 



(8) 



where we may term /(r) the 'profile' while the phase (p 
coincides with the azimuthal angle and n is an integer. 
Using the GP equation ^ yields 



/ 



ld_ 

r dr 



dr 



f + I\f\\ 



(9) 



which can be integrated numerically in an infinite 
domainii^ to determine the nature of f(r). The bound- 
ary conditions require that /(r) — > 1 as r — oo, while, in 
the core of the vortex, f{r) as r 0. 

Since, as mentioned, there is no analytical expression 
for /(r), various approximate formulae have been pro- 
posed. Among these are (i) a [1/1] Padc approximant 
advanced by Fetter;^ namely. 



fF{r) 



(10) 



with p = 2; (zi) a 'square core' approximant, introduced 
by Kerr,— say /k, also given by (|T0| but with p = 4. 
This particular choice was introduced as a stratagem to 
ensure that the initial density went smoothly from zero, 
on the vortex core, to roughly the background density 
over the distance of the healing length Finally, (iii) 
BcrlofS^ derived a [2/2] Padc approximant specified by 



/B(r-) 



r2(0.3437-f 0.0286r2) 
1 -I- 0.3333r2 -|- 0.0286r4 



(11) 



A precise stationary solution, which for convenience we 
will call the diffused profile fair), can be found by evolv- 
ing in time through the DGP equation any one of these 
three profiles. Since we are imposing periodic boundary 



conditions, we need a large enough domain to minimize 
the infiuence of the images. We have calculated /^(r) 
in a domain of size ~ Ly = 200 using a grid spacing 
Ax = Ay ~ 0.098. 

To improve on the options we have generated a two- 
point [2/2] Padc approximani»^ in the form 



/[2/2] ('') 



ciQ + air- + 



[2/21^'^ 1 + 5ir2 -I- 62r4 
and a corresponding [3/3] Padc approximant 



/[3/3](^) 



Co + cir^ + C2r^ + c^^r^ 



1 



(12) 



(13) 



+ dir"^ + d2r^ + dy.r^ 
By recalling that, at small r, the profile /(r) can 
be expanded as /(r) ~ X^i^i P^'^^*^^' with^i pi ~ 
0.5827811878, P2 = -pi/8, and = Pi(p? + l/8)/24, 
while about ?- — > oo one has^i /(r) ~ Y^^n QiT'^'^ . with 
go = 1, gi = -1/2 52 = -9/8 and 93 = 
find, for /[2/2], Oo = 0, oi = p?, 02 = 62 



= -161/16, we 
and 



hi = 



4(1 

^2 =pI 



-vl) 



0.3857, 



4(1 



0.0461. 



(14) 



The expressions for hi and 62 are found by expanding (jl2p 
in a Taylor series for small r and large r and equating the 
coefficients of the second term, respectively, to p2 f^i^d 
qi. This approximant then reproduces pi, P2 md qq, qi 
correctly. 

Following an analogous procedure we find, for /[3/3], 
Co = 0, ci = pI, C3 = ds, C2 ~ 0.0501, di ~ 0.3976, 
d2 — 0.0527 and ^3 ~ 0.0026. This approximant repro- 
duces pi, P2, P3 and go, 91, 92- The /[2/2] and /[3/3] Fade 
approximants and the Berloff approximant are compared 
in Fig. [T] with respect to the diffused profile fo- The 
[3 /3] Pade approximant is the closest to the diffused pro- 
file followed by the [2/2] Pade approximant and, lastly, 
the Berloff approximant. 

It is observed that expanding fp, fx and Jb in power 
series, and examining the limits for small and large r, 
yields a variety of results, none of which reproduces pi, 
P2 and qi. Explicitly, for small r, one has 



h 



r 

^2 



4V2 

^6 



+ 0(r^), 



^ — - — + 0(r^-') 



^/2 4V2 



(15) 



V0.3437r - 0.2501y VO.3437 + 0{r^); 



while the behavior for large r is 



1 



If 

fK 



1 



O 



llj.2 +^(7.4 



(16) 
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FIG. 1. Difference between the Berloff approximant and the diffused solution (dash-dotted line), between the [2/2] Fade 
approximant and the diffused profile (dotted line), and between the [3/3] Fade approximant and the diffused profile (dashed 
line). Inset: time decay of the ratio between the total energy E and the total initial energy Eq for the [2/2] Fade approximant 
evolved through the DGP equation. 



Profile 


Lo 


fF 


0.403 


fK 


0.604 


fB 


0.390 


/p[2/2] 


0.382 


/P[3/31 


0.381 


fD 


0.377 



TABLE I. Lo parameter for the five different approximate profiles and the numerical diffused profile. 



Notice that, even though /s — > 1 as ?- — > oo, the BerlofF 
profile, /s, exceeds unity somewhat in an intermediate 
range of r, which is qualitatively incorrect. 

By imposing the Fetter and Kerr profiles as initial con- 
ditions and by evolving them through the GP equation 
in a periodic domain it is also discovered that, as a con- 
sequence of being less accurate approximants, they are 
markedly unstable: see Fig. 2 of Meichle et al,^ In- 
deed, the core region first relaxes toward the exact pro- 
file while the excess energy released is successively prop- 
agated outwards as waves. These remain localized in the 
outer region of the domain, where they persist and are 
responsible for background oscillations which one might 
be tempted to identify as 'thermal noise'. If these vortex 



profiles are imposed for two- vortex initial conditions, and 
vortex annihilation or vortex reconnection are studied, it 
is observed that the presence of these background 'acous- 
tic' waves strongly affects the vortex trajectories and the 
speed of the processi^ In particular, vortex trajectories 
become more wavy and annihilation or reconnection take 
longer to occur. One may speculate that in the presence 
of a 'thermal background' some effective frictional effects 
arc generated; wc comment on this further in the follow- 
ing section. 

Being a stationary solution of the GP equation, the ex- 
act single vortex profile minimizes the Hamiltonian. For 
this reason the value of the total energy is indicative of 
the accuracy of the approximate profiles. The energy per 
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unit length of a single vortex with n = 1 in a cylindrical 
domain of radius R is given hy^ 



E' =TT- 





\df{r)^ 


{r 


dr 



rdr 



fir) 



dr 



(17) 



where the first term gives the quantum energy E'q, the 
second, the classical kinetic energy E'^^^^ and the third, 
the potential energy E'j. 

For large R we have /(r) ~ 1 and P7)) simply gives 
E' sa Trlni?. For greater asymptotic precision Pitaevskii^ 
found E' = 7r[ln(i?) + Lq], with Lq = 0.38. The values 
of the parameter Lq for all the five profiles and for the 
numerically diffused profile are reported in Table HI 



IV. VORTEX DIPOLE DYNAMICS 

Customarily, as mentioned in the Introduction, initial 
conditions for multiple vortices have been constructed 
simply by taking the product of the wave functions for 
each individual vortex Ji^ii that is 



vl/o(x,t = 0) = n/je 



(18) 



where fj (x) and (pj (x) are the profile and phase of the 
j-th vortex with specified core-locus. In the case of a 
vortex dipole characterized by a separation distance do, 
diff'erent two-vortex profiles provide different choices of 
fj [see Fig. [2](a)], while the sum of the phases <f>j leads 
for all cases to the 'circular' phase contours, illustrated 
in Fig. [2jb). We have, in addition, taken account of the 
periodic images in accord to the explanations at the end 
of Sec. II above.— 

In our vortex dipole calculations we have used a 2D 
domain of size 
range 



Ly ~ 50, spanning the coordinate 



25 ^ a; > 25 , - 25 ^ y ^ 25, 



(19) 



with X ~ and \y\ = a ~ do/2 setting the location of 
the vortex cores. The initial global phase (po, according 
to (UHl, is then 



(f>o{x,y;t = 0) = 01 + </>2 

1 / w — a 
^ tan~^ ' ^ 



tan 



y + a 



tan 



-lax 



x^ ^ — 



(20) 



Now consider an alternative approach in which: (i) 
the phase contours (PD|) are generalized to the elliptical 
distributions 



4>o{x, y;t = 0) ^ tan 



-2aux 



K[{ux)^ + y'^ — a? 



(21) 



with u and K as parameters; (iz) the profiles for a trial 
initial condition are obtained using (jl8p , but are relaxed 
through the application of (l7|), namely the DGP equation 
[see Fig. [21 parts (c) and (d)], while the phase is fixed by 
([2T|). Clearly the values K = u = 1 correspond to the 
'circular' profile given by (PU)) . 

In (PT|) the loci of fixed 0o are ellipses that pass through 
the dipolar cores at (x, y)—{Q, ±^do) with axes parallel 
to the x and y axes. When (po = ±7r/2 the ellipse is 
centered at the origin and has width xo = ±(a/u). If u — >■ 
this ellipse degenerates into a pair of parallel horizontal 
straight lines. When 0o is small the center of the ellipse 
is displaced to xo ~ —l/Kucpo and the width increases 
by a factor -^/l + I/^qA'w. Finally, for u < 1 the ellipse 
is elongated along the x axis, but for u > \ along the y 
axis, while, when K < 1 the phase contours at fixed (po 
are further apart than when A' = 1. See Fig. [3] for an 
illustration. 

Note that two antiparallel vortices forming a vortex 
dipole are not a stationary solution of the CP equation. 
Hence if such a configuration is allowed to diffuse fully 
through the DGP equation, it converges to the fiat solu- 
tion, 4* = 1. However, when phase contours are enforced, 
i.e., held fixed while diffusing, a two- vortex combined 
profile relaxes to a lower energy solution which preserves 
the presence of both vortices as illustrated in Fig. [5J 

In the split-step spectral method each time increment 
At is divided into two steps: the first integrates the non- 
linear term in the Cartesian frame of reference, the sec- 
ond integrates the linear term in spectral space. The de- 
sired phase 00 is imposed by redefining, after each step, a 
wave function ^enforced phase given by the modulus, |*|, 
of the wave function just calculated times the exponential 
of 00 as specified in ((2T|) so that VPenforced phase = l^le*"^". 

A hint as to the outcome one may expect from the time 
evolution of such initial conditions is obtained by analyz- 
ing the linearized CP equatiouj^i which, in fact, is sim- 
ply the time-dependent one-particle Schrodinger (OPS) 
equation. By applying a phase shift ^' = $e**/^ and ne- 
glecting the nonlinear term, the CP equation reduces to 
the OPS equation which, in terms of <l?(x, t), is 



2i 



'dt 



= 0. 



(22) 



Given the elliptical phases (f2T|) the chosen initial con- 
dition becomes 

{uxf + - + 2iaxu/K 
*o T^) ^^^^ 

where V{r) is any appropriate polynomial that allows for 
the correct asymptotic behavior of 1^*0 1 for r close to and 
far from the vortex cores. For the sake of illustration 
and to simplify the calculations we take T' = 1; in doing 
so, however, the correct asymptotic behavior is not guar- 
anteed, especially for r far from the vortex cores since 
l^*! — >■ oo as r — > oo. In (|23|1 the imaginary and real 
parts of 5*0 are simultaneously zero at (x = 0, y — ±a), 
which, indeed, provides the location of the vortex lines. 
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FIG. 2. (a) Two- vortex profiles viewed along (x = 0, y), for do = 2.5, obtained by multiplying wave-functions for each of the 
two vortices as in (|18p : the Fetter profile (green, dashed line), the Kerr profile (black stars) the Berlofi' profile (blue circles); 
or, by relaxing through the DGP equation any one of the three previous two-vortex profiles while keeping the initial phase 
contours fixed (red, solid line); (b) phase contours given by (|20p . associated to the profiles shown in (a) and plotted every 7r/4; 
(c) total energy per unit mass for the two- vortex Fetter profile (diamond), Kerr profile (star), Berloff profile (circle), and time 
decay of the total energy observed letting the DGP equation evolve from the two- vortex Fetter profile (dotted line). In this 
plot the maximum value is normalized to unity; (d) density distribution, i.e., |^'|^ contours, for the diffused profile. 



The solution of (g^l) with $(x, <=0) = *o, is found to 
have real and imaginary parts given by 



Re($) = (uxf 



2 2 

V -a , 



Im($) = 2axu/K + t{u^ 
which simultaneously vanish when 

X = -t{u'^ + l)K/2au, 



1), 



2/1,2 = ±1 



t(u2 + 1)K 



2a 



(24) 



(25) 



Hence, in this limit, the vortex cores move linearly with 
time in the direction of negative x along the ellipse which 
coincides with Re(^'o) = 0. On equating the first space- 
derivative of |<f>(x, t)p to zero one finds two solutions: 



the first reproduces and corresponds to the moving 
cores, i.e., the zero-density points; but the second one is 
given by y = and 



u'^x^ 



+ a^u^x{2/K - 1) + autiu"- + l)/K = 0, (26) 



which can be solved using Cardano's formula for the roots 
of a cubic equation. The real root of this equation cor- 
responds to the location of the minimum density points 
after annihilation of the vortex dipolc. 

We conclude that, given the elliptical phase initial con- 
figuration, a vortex dipole with V constant or of first 
order in {x,y) should be expected to annihilate via the 
linearized GP equation in accord with (P^. If V is of 
higher order, (j24p represents an approximate solution ac- 
curate to first-order in time. The approach that lead us 
to find solution (|24p has been proposed and followed by 
Nazarenko and West<^ 
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FIG. 3. Phase contours given by eq. H2ip and plotted every 7r/4 for difTerent values of the eUiptical parameters u and K. 



Observe that if m — >■ and uw = K ^ 0, that is, u 
and K go to zero, but in fixed proportion with w as the 
control parameter, the condition ([23|) becomes 

*o = (y^ -a^) +2ma;w;. (27) 

The OPS solution for this is 

$(x,t) = - + 2aa;w). (28) 

The time behavior of the vortex core is then given by 

2/1.2 = ±a and x = —t/2aw. (29) 

Thus the vortex dipole now does not annihilate but, 
rather, moves with a fixed velocity of propagation in- 
versely proportional to a and w. 

In the following we solve numerically the nonlinear GP 
equation for different vortex-dipole initial conditions and 
we show: (a) how annihilation of a vortex dipole may 
be affected by the choice of different initial profiles fj\ 
(b) how, given 'circular' phase contours and a diffused 
two-vortex profiles, the solution depends on the initial 



separation distance do] and (c) how, the vortex-dipole 
dynamics is influenced by varying the initial phase con- 
figurations via the elliptical parameters u and K in (|2ip . 



A. Sensitivity to choice of initial profiles 

In the first case we take K = u — I in (PT|). i.e., cir- 
cular phase contours, and examine different two- vortex 
initial profiles. We find that if the diffused structure (at 
fixed phase distribution) for a vortex dipole is imposed 
as an initial condition for the GP equation, the annihi- 
lation mechanism proceeds at a markedly different rate 
than when the initial condition is obtained by multiplying 
the two single- vortex profiles: see the open-circles plot in 
Fig. m Evidently the diffused profile commences annihi- 
lating sooner and the annihilation process proceeds more 
rapidly than do the other three profiles. This may be 
rationalized by noting that the non-diffused profiles pos- 
sess excess energy that needs to be radiated away from 
the cores before the annihilation itself can get properly 
underway. 
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FIG. 4. Time evolution of the minimum density for two an- 
tiparallel vortices with initial separation do = 2.5, starting 
from four different initial conditions: the two-vortex Fetter 
(violet dots), Kerr (black stars), Berloff (blue crosses), and 
[3/3] Fade (green triangles) profiles and any of the previous 
conditions diffused through the DGF equation, (0, keeping 
the phase contours fixed (red open circles). 
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FIG. 5. Position in the (x, y) plane of the vortex cores as 
time proceeds, comparing the cases when the initial profiles 
are derived from the Kerr (stars), Berloff (crosses). Fetter 
(dots), [3/3] Fade approximants (diamonds), with the relaxed 
or diffused two- vortex profile (circles). 



In confirmation of this interpretation, the Kerr pro- 
file is the one which starts the annihilation process most 
slowly, followed in order by the Fetter, Berloff, [3/3] Fade, 
and diffused profiles. The detailed trajectories of the 
vortex centers before annihilation are also strongly influ- 
enced by the choice of the initial condition as illustrated 



in Fig. [5] for the case do = 3. For these simulations the 
conservation of the total energy AE is accurate to within 
0.3%. 

The dipole annihilation mechanism can be thought to 
proceed generally in two stages: (i) the two vortices 
translate and move towards one another while the phase 
reduces its range of variation everywhere except close to 
the vortices: see Fig. [Sfa) and (b); {ii) the cores merge 
together, annihilation occurs and the minimum density 
increases while the whole structure continues to move 
forward radiating acoustic waves: Fig. ^c) and (d). 
This description is consistent with the results reported 
by Ogawa et al^ 

By choosing a two-vortex diffused profile and repeat- 
ing the calculations for different initial separations, it is 
found that annihilation always occurs when c?o ^ 3. If, 
instead, the initial separation exceeds dp ~ 4 the vor- 
tex dipole eventually propagates as a solitary wave (al- 
beit usually with associated oscillations). It transpires, 
however, that this observation is deeply connected with 
the specific choice of the initial phase distribution. We 
note, indeed, that the phase contours entailed in the ap- 
proximate solitary wave solutions found by Berlof]^ for 
do = 1.8 are more-or-less elliptical in shape, rather than 
circular, as those imposed here: see Fig. Elja). Inter- 
estingly, when our do = 4 initial condition is evolved 
through the GP equation, it behaves as a solitary wave 
and its phase distribution, which is initially circular as in 
Fig. ^c), relaxes to an elliptical- like configuration rather 
similar to Berloff 's result as shown in Fig. [7{b). 



B. Dependence on the initial separation distance do 

Starting from circular phase contours {K = u = 1) 
and a two-vortex diffused profile we first report in Fig. 
[5I^a) the the annihilation time, tji^a{do), needed to just 
complete the merging of the two distinct cores. The ap- 
proximately quadratic dependence of iAn on do observed 
might have been anticipated from the solution ([25]) of the 
linearized GP (or OPS) equation (22). Indeed, this yields 
an annihilation time 

tl^{do) = bodl with bo = l/2Kiu^ + l), (30) 

where the superscript denotes the linearized derivation, 
despite which the result seems likely to be generally valid 
for do 0. Accordingly in Fig. ^a), where K = u = 1, 
we present a fit using 



*An(c^o) - ^d,l[l + bido] 



(31) 



with 6i = 0.39. Clearly these limited data are consistent 
with the form pop : however, while recognizing the rela- 
tively large uncertainties involved in these calculations, 
it should, perhaps, be mentioned that a direct estimate 
of bo yields a somewhat larger value around 1/3. 

Once the vortex dipole has been annihilated there re- 
mains only a 'dip' characterized by a minimal density. 
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FIG. 6. Annihilation process through the GP equation for do = 1 showing only a small portion of the domain. At t — 0.25: (a) 
distribution of (b) phase contours at multiples of 7r/20. At t = 2: (c) distribution of |>1/|^ after annihilation, the minimum 
density exceeding zero; (d) phase distribution now confined to the interval —n/20 < < pi/5. 



I^(*)lmin' ''^hich, as t grows above increases mono- 
tonically from zero but initially fairly slowly: see the 
plots in Fig. [Hb), especially those for dg > 1-5. Then, for 
larger times, the density minimum evidently approaches 
a plateau. The process depends strongly on the initial 
separation do in a manner that seems linked to the total 
missing mass associated with original vortices, namely, 
/^(l — Furthermore, the visually quite varied 

plots in Fig. [H^b) can be reduced to a single 'master 
growth curve', as seen in Fig. |Hl[c), by rescaling the 
time by tAn('^o) and the density by a characteristic den- 
sity I ^' (do) I To determine I^'Iaq nsed the value 
\^{do;t = 2iAn)P which, as reported in Fig. [8l[d), ini- 
tially rises linearly with do with a slope (indicated by a 
dashed line) of about 0.033. 



C. Sensitivity to initial phase configuration 

As our last study a set of initial conditions, with pro- 
files relaxed via the DGP equation, are imposed by en- 



forcing the elliptical phase contours stated in (|2ip for a 
range of the parameters u and K. The resulting transi- 
tion between the annihilating and propagating behavior 
of the vortex dipoles is evident in Fig. |9] for the three 
separation values: (a) do = 1.8, which is the separation 
for which Berloff derived a solitary- wave approximantj^ 
(b) do = 2.7 and (c) do = 4.0, the value for which we 
originally observed solitary-wave behavior with initially 
circular phase contours. The annihilation-solitary wave 
(or AnSol) boundary clearly depends on both the ellip- 
ticity parameters, K and u, and on the separation do 
between the parallel vortices. 

Note, first, that the data suggest quite strongly that 
the AnSol boundary approaches the {K, u) origin linearly 
when u (and K) become small. In this limit we find the 
boundary is reasonably well represented by 

Kc(u;do) ~ cowdp, (32) 

with Co ~ 0.056 for do < 4. On the other hand, a 
theoretical justification of this expression even, say, for 
do — !• (for example on the basis of solutions of the lin- 
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FIG. 7. (a) Phase distribution of the solitary-wave approximant derived by Berlofl^ for steady velocity 11— OA, to be compared 
with (b) the phase distribution at t=20 for an initial condition characterized by do = 4 after evolution by the GP equation. 




FIG. 8. Vortex dipole GP annihilation dynamics for circular phase contours {K — u — 1) relaxed via the DGP equation, (a) 
The annihilation time tj^j^{do), and a simple fit, quadratic as do — >■ 0: for details see text, (b) The minimum density Ivl'lmin 
a function of time as the vortex dipoles merge; (c) the collapse of the |^|min plots rescaled by tAn('^o) and by the characteristic 
density, shown in (d) as a function of do. See further in the text. 
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FIG. 9. Annihilating and propagating behavior for a vortex dipole characterized by elliptical phase contours specified by {K, 
u) for separations: (a) do = 1-8, (b) do = 2.7, and (c) do = 4.0. Dots (•) correspond to annihilating dipoles, and crosses (+) to 
propagating dipoles. The circles (o) indicate calculations for parameters close to the AnSol boundary in which annihilation is 
almost observed, as the vortex cores approach closely, although, subsequently, a strongly oscillating propagating dipole appears 
to develop. The solid lines (color online) are guides to the eye as regards the AnSol boundaries. 
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earizcd GP equation or the corresponding solitary wave 
GP equation3^) appears nontrivial. 

Then on noting the distinct vertical scales in Fig. [9] 
we see that for large u the AnSol boundary becomes in- 
dependent of u so that Kc(u; do) — A'oo(rfo) as u — > oo. 
However the larger the vortex separation do, the larger 
the value of K at which the AnSol transition is realized. 
In fact it seems that 

i^oo(do) — Coodo, (33) 

with Coo — 0.091 provides, within rather large uncertain- 
ties, a tolerable description. 

Lastly, it appears that at least for do > 2, the AnSol 
boundary is nonmonotonic in u; rather Kc{u; do) displays 
a maximum around u = 1-2 which, relative to A'oo(do), 
also increases quite dramatically with do. More extensive 
calculations would be needed, however, to provide more 
quantitative conclusions. 

V. CONCLUSIONS 

A single straight vortex is usually regarded as a station- 
ary axisymmetric solution of the Gross-Pitaevskii equa- 
tion. We have demonstrated, however, that some approx- 
imate vortex profiles proposed in the literature contain 
a significant excess of energy which, under time develop- 
ment, leads to the outwards propagation of compressional 
waves. These effects prove quite small for the particu- 
lar Pade approximant proposed by BerlofS^ and are cer- 
tainly negligible for our newly proposed exact two-point 
[2/2] and [3/3] Pade approximantS)^ Thus these profiles 
qualify as good candidates to perform low-energy calcula- 
tions: a desirable condition, since the background 'ther- 
malization' which arises from less accurate profiles is an 
uncontrolled and, initially at least, an unwanted effect. 
Furthermore, the undetermined excess energy interferes 
with the precise quantitative measurements one should 
hope to perform. It may be expected, moreover, that 
real physical systems in which energy dissipation (not 
included in the GP model) is to be anticipated, will tend 
to assume minimal energy configurations when feasible. 

The second part of our study has been devoted to the 
dynamics of vortex dipoles, characterized by a separation 
distance do between two parallel counter-rotating infinite 
straight line vortices. It was demonstrated that such a 
dipole can either annihilate or propagate steadily as a 
solitary wave. Both the qualitative and the quantitative 
behavior depend strongly on the initial vortex profile and 
on the initial phase contours. Starting from an initial 
condition obtained through multiplication of two single- 
vortex wave functions one typically observes annihilation 
for an initial separation, do, smaller than three intrinsic 
units, but solitary-wave behavior for an initial separa- 
tion larger than four units. When do is small, a good 
approximation for the dynamical behavior is given by 



the linearized GP equation. The linear term then dom- 
inates and annihilation is observed, as predicted by the 
one-particle Schrodinger equation. If lower energy ini- 
tial profiles are used, annihilation occurs faster. The do 
dependence of the solutions has been explored by deter- 
mining the low-energy initial conditions via the Diffuse 
GP equation and the time evolution via the GP equation. 

When the separation distance is larger, the linearized 
solution is no longer valid, the nonlinear term plays a role 
and its interplay with the linear terms determines the 
solitary wave behavior; of course, this is a typical non- 
linear phenomenon as observed, e.g., in the Korteweg-de 
Vries equation . ^^'^^ It is interesting to note that, in this 
case, starting with a rough, approximate solitary-wave 
as the initial condition, the solution rapidly descends, 
as time develops, into the 'solitary-wave valley' in an 
idealized landscape map of possible solutions. On the 
other hand, when the initial separation distance is small, 
this outcome is rarely observed. However, solitary wave 
solutions for small distances, even for the limiting case 
do = 0, do exist; but they require precise initial condi- 
tions such as those found by Berloff<^ 

More quantitatively, we have thus studied a class of 
explicit initial phase configurations that yield elliptical 
phase contours specified by two parameters, K and u: see 
the expression (PT|) and Fig. [5] which reveals well-defined 
annihilation vs. solitary wave (or AnSol) boundaries of 
characteristic form. 

In conclusion, we remark, in more general terms, that 
vortices are specified by the intersections of the zeroes of 
the real and imaginary parts of the effective wave func- 
tion: these form 'zero lines' (or linear loci) in two spa- 
tial dimensions, or 'zero planes' (or planar manifolds) in 
three dimensions. There are many possible topologies, 
and functional forms associated with them, that allow 
for vortex intersections. Fixed points of the GP equa- 
tion then correspond to well defined topologies;^ but for 
dynamical structures such as a vortex dipole or recon- 
necting vortices the choice of initial configurations re- 
mains arbitrary and may strongly infiuence the outcome. 
One must conclude that the custom of representing the 
various topologies merely by the product of approximate 
wave functions derived from straight vortices is likely to 
be a significant limitation in developing a valid perspec- 
tive for understanding the entire dynamic landscape more 
fully. 
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